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ABSTRACT 



Context. Far-infrared dust emission has a self-similar structure which reveals the complex dynamical processes that 
shape the interstellar medium. The description of the statistical properties of this emission gives important constraints 
on the physics of the interstellar medium but it is also a useful way to estimate the contamination of diffuse interstellar 
emission in the cases where it is considered a nuisance. 

Aims. The main goals of this analysis of the power spectrum and non- Gaussian properties of far-infrared dust emission 
are 1) to estimate the power spectrum of interstellar matter density in three dimensions, 2) to review and extend previ- 
ous estimates of the cirrus noise due to dust emission and 3) to produce simulated dust emission maps that reproduce 
the observed statistical properties. 

Methods. To estimate the statistical properties of dust emission we analyzed the power spectrum and wavelet decom- 
position of 100 /im IRIS data (an improved version of the IRAS data) over 55 % of the sky. The simulation of realistic 
infrared emission maps is based on modified Gaussian random fields. 

Results. The main results are the following. 1) The cirrus noise level as a function of brightness has been previously 
overestimated. It is found to be proportional to < / > instead of < / > 15 , where < / > is the local average brightness 
at 100 /im. This scaling is in accordance with the fact that the brightness fluctuation level observed at a given angular 
scale on the sky is the sum of fluctuations of increasing amplitude with distance on the line of sight. 2) The spectral 
index of dust emission at scales between 5 arcmin and 12.5° is < 7 >= —2.9 on average but shows significant variations 
over the sky. Bright regions have systematically steeper power spectra than diffuse regions. 3) The skewness and kurtosis 
of brightness fluctuations is high, indicative of strong non-Gaussianity. Unlike the standard deviation of the fluctua- 
tions, the skewness and kurtosis do not depend significantly on brightness, except in bright regions (> 10 MJy sr _1 ) 
where they are systematically higher, probably due to contrasted structures related to star formation activity. 4) Based 
on our characterization of the 100 /im power spectrum we provide a prescription of the cirrus confusion noise as a 
function of wavelength and scale. 5) Finally we present a method based on a modification of Gaussian random fields 
to produce simulations of dust maps which reproduce the power spectrum and non- Gaussian properties of interstellar 
dust emission. 

Key words. Methods: statistical - ISM: structure - Infrared: ISM: continuum - dust 



1. Introduction 

The interstellar medium emission shows fluctuations at all 
observable scales, revealing the self-similar nature of its 
density structure. The physical processes responsible for 
this self-similarity of the interstellar medium (ISM) struc- 
ture is yet to be fully identified. It could be related to turbu- 
lent motions but also to chemical and thermal instabilities 
which trigger phase transitions and play an important role 
in shaping the medium. 

Faced with the challenge to understand and char- 
acterize the great structural complexity of interstellar 
emission, astronomers have used several statistical tools 
(power spectrum, correlation and structure functions, 
wavelets, area-perimeter relation, principal component 
analysis,...) on several interstellar tracers : mo l ecules 
( Falg arone et all 119911: iHobsonl 1992; Stutz ki et ali i 1998: 



Stanimi rovic et all 

~| Dickev et all 



119991: 
2001; 



Stanimirovic & Lazarian, 



2001; Dic kev et all 120011: lElmegreen et al 
Miville-Desch enes et alll2003af) . extinction (iPadoan et al 
20021 " 



1992 



2004). 



Pad oan et al .. 2006) or d ust em i ssion (IGautier et al 



. _ , „„„ __ {& 

lAbergel et all . Il996l : LlewellL l200ll : Engalls et al 



iBensch et a ll l200lt iBruntl 120031: IPadoan et all 12003) , 
atomic hydrogen ( Crovisier fc Dickevl Il983l : lGreenU l993; 
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The analysis of the results of these tools is made 
difficult by projection, and instrumental effects but also 
by the fact that no observation is a perfect tracer 
of the total gas column density. Several works have 
been dedi c ated t o the understan d ing o f these effects 
(iGoldmanl. l2000l: lElmegreen et all |200ll : IPadoan et all . 
120011 : iMiville-Deschenes et all l20Q3bh and specifically on 
how one can retrieve the three-dimensional statistical prop- 
erties of interstellar matter from astronomical observations. 
The recent theoretical progresses made in that area open 
interesting perspectives on the use of statistical tools to 
determine the three-dimensional structure of the gas and 
to identify privileged scales at which physical processes are 
important. 
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l 10 „ (MJy/sr) 

Fig. 1. Top. Cumulative histogram of 100 /im surface 
brightness over 98% of the sky (solid line): fraction of the 
sky with a brightness lower or equal to iioo- The dashed 
line represent what would be expected from a cylindrical 
disk (cosecant law). Bottom. Histogram of iioo bright- 
ness over the whole sky with bin scaled logarithmically. For 
both curves we used the IRIS data projected on the Healpix 
grid (equal area pixels - http://healpix.j pT.nasa.gov] ) with 
a pixel size of 1.7 arcmin (nside=2048). A background of 
0.78 MJy sr -1 was removed from the data to take into 
account the co smic infrared background and any zodia cal 
light residual (|Hauser et all Il998cl lLagache et all [2000V 



The characterization of the statistical properties of the 
interstellar emission is relevant for the understanding of 
the physics of the ISM but is also of importance for studies 
where interstellar emission is a nuisance. This is the case of 
the study of pre-stellar cores in molecular clouds but also of 
the analysis of the Cosmic Microwave Background (CMB) 
and Cosmic Infrared Background (CIB) radiations. In these 
cases the interstellar emission is considered as a noise (the 
so-called "cirrus noise" ) for which one needs a detailed sta- 
tistical description in order to remove it and account for it 
in the error budget. In that respect the interstellar medium 
is a rather complex noise source due to its highly non- white 
and non- Gaussian brightness fluctuations. 



< iioo > 


16.18 MJy sr~ 


median iioo 


3.49 MJy sr~ 


most probable iioo 


2.51 MJy sr" 


< 0.1 MJy sr -1 


0.006% 


< 0.2 MJy sr -1 


0.06% 


< 0.5 MJy sr -1 


1.8% 


< 1.0 MJy sr -1 


11.6% 


< 2.0 MJy sr -1 


29.5% 


< 5.0 MJy sr -1 


61.1% 


< 10.0 MJy sr -1 


77.6% 


< 20.0 MJy sr -1 


89.2% 


< 50.0 MJy sr" 1 


95.6% 



Table 1. Average, median and most probable value of 
the 100 /im brightness over 98 % of the sky at an angular 
resolution of 4.3 arcminutes. The bottom portion of the 
table gives the fraction of the sky with brightness lower 
than the value given in the left column. 



The present paper is a study of the power spectrum of 
interstellar dust emission at 100 /im. The main goals are 
to use such analysis to put some constrains on the den- 
sity structure of the interstellar medium but also to pro- 
pose a caveat f or est imation and simulation of cirrus noise. 
iGautier et al.l ([1992) showed that the power spectrum of 
the IRAS 100 /im emission is characterized by a power law 
P(k) = Ak 1 where the e xponent 7 ~ — 3 is in dependant of 
the brightness. Moreover IGautier et al.l (|l992l ) showed that 
the normalisation factor of the power spectrum depends 
on the mean brightness < I > of the region considered, 
with A oc< I > 3 . This study has been conducted on a rel- 
atively small number of regions and at a time where the 
instrumental response of IRAS was not well known. In this 
paper we would like to revisit the work of IGautier et al.l 
(1992) by extending it to the whole sky and by taking ad- 
vantage of the recent reprocessing of the ISSA plates by 
IMiville-Deschenes fc Lagachd (|2005f ). 

In §[2] we present the data used in this analysis. The re- 
sults of the power spectrum analysis and the implication on 
the ISM density structure are presented in § [3] An estimate 
of the cirrus noise level in the FIR-submm is provided in 
§ 3] and a method to produce realistic dust emission maps 
is given in § [5j 

2. The IRAS survey and IRIS 

In 1983 the Infrared Astronomical Satellite (IRAS) has 
made a survey of 98% of the sky in four bands : 12, 25, 
60 and 100 jam. Since then this dataset has been used ex- 
tensively in almost all area of astrophysics. In the early 
90s the IRAS team released the IRAS Sky Survey Atlas, 
a set of 430 fields with better calibration. Each field is a 
12.5° x 12.5° image with a pixel size of 1.5' x 1.5'. Recently 
IMiville-Deschenes fc Lagachd (|2005f ) reprocessed the ISSA 
maps to correct for residual calibration defects and strip- 
ping. This new set of ISSA plates, named IRIS, improves 
significantly the quality of the IRAS data by lowering the 
noise level and improving the photometry of the four bands. 
Our analysis is based on the 100 /im IRIS maps, available 
in native 12.5° x 12.5° cartesian maps or in H ealpix vector 
at |http: / / www.cita.utoronto.c a/ ^mamd/IRIS / . 

One important parameter of our analysis is the esti- 
mate of the instrumental noise contribution to the power 
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Fig. 2. Power spectrum of a typical IRIS map (point 
sources removed) with its associated noise power spectrum 
and the estimated CIB level (convolved by the IRAS beam). 
The power law at large scale (small k) is due to the Galactic 
dust emission. The cutoff at k ~ 0.1 arcmin -1 is due to the 
IRAS beam. 




1 10 100 

<l> (MJy/sr) 



Fig. 3. Spectral index 7 of the power spectrum for our sam- 
ple. Top. Histogram of the spectral index values. Median 
is -2.93, average is -2.96 and standard deviation is 0.3. 
Bottom. Spectral index as a function of average 100 jam 
brightness (CIB and zodiacal light subtracted). 



spectrum. We used only the IRIS maps for which we could 
compute the contribution of the noise. To do so we used 
the fact that the original ISSA plates are the combination 
of up to three HCON images. The power spectrum of the 
noise can be estimated by taking the power spectrum of 
the difference b etween two HCON images, as described in 
more details by IMiville-Deschenes et all ((2002). 

We selected only maps for which each pixel has been 
observed at least two times so that a noise map can be com- 
puted, which represent 236 maps (out of 430). Most maps 
of our sample have an average brightness < I > lower than 
20 MJy sr -1 . This is representative of the whole sky statis- 
tics, as seen in Fig. [1] (see also Table [1]) where the proba- 
bility density function (PDF) and cumulative histogram of 
the 100 jam brightness (CIB and zodiacal light subtracted) 
over the whole sky is presented. This figure, together with 
Table [TJ show that 90 % of the sky has a brightness lower 
or equal to 20 MJy sr -1 (or Nh ^ 4 x 10 21 cm -2 according 
to lLagache et al.l (|2QQ0h ). This simple statistics of 100 jam 
brightness reveals that there is less than 2 % of the sky 
with N H < 1 x 10 20 cm" 2 . 



3. Power spectrum analysis 

The power spectrum has been extensively used in the anal- 
ysis of image structure. Working in Fourier (or spherical 
harmonics) space offers considerable advantages; this for- 
malism allows for the quadratic separation of uncorrelated 
sources in the image (e.g., noise and signal) and decon- 
volution of the instrumental function. In addition working 
in Fourier space facilitates the comparison with numerical 
simulations of turbulent flows. 



3.1. Computation of the power spectrum 

The power spectrum of an image f(x,y) of Fourier 
Transform f(k x ,k y ) is computed from the amplitude 
A(k x ,k y ) defined as 

A(k x , ky) = / (k x , ky) f (k x , k y ) = \f(k x , k y )\ . (1) 

The power spectrum P(k)dk is an angular average of 

A(k x ,k y ) between k and k + dk where k = yjk%. + ky. The 

method we u se to compute the power spect rum is the one 
described by IMiville-Deschenes et all ([2002). To minimize 
edges effects in Fourier space we apodize the image, from 
which the mean was removed, using a cosine tapper of 15 
pixels wide. 

3.2. Contributions to the power spectrum 

The main goal of our analysis is to characterize the power 
spectrum of the interstellar diffuse emission, but several as- 
trophysical signals and artifacts may contribe to the power 
spectrum of the IRIS maps. The power spectrum P{k) of 
the IRIS maps may be formalized by the following equation: 

P(k) = B(k)(P ism (k) + 

Psources (k) + P cib ) + N(k) (2) 

where B(k) is the effective beam of the IRIS maps (which 
includes the instrumental beam and the map making), 
P%sm(k), Psources (k), P cib (k) and N(k) are respectively the 
contributions of the interstellar medium, detected point 
sources, the unresolved CIB and the noise. All these contri- 
butions have to be estimated to characterize the interstellar 
component. 

A typical power spectrum of a faint region of the sky is 
shown in Fig. [21 At large scales (low k), the power spectrum 
follows a power la w Ak 1 ', this is the si gnature of the cir- 
rus cloud emission (jGautier et al.l .[l992). At smaller scales, 
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Fig. 4. Normalisation P(0.01) of the power spectrum at 
0.01 arcmin -1 of each map as a function of its mean 
100 /im brightness. The solid line is our fit to the data 
(see Equa tions [5] and 6 t. Th e dashed line is the relation 
given bv lGautier et al.P|1992[ ). 
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Fig. 5. Standard deviation as a function of the average in- 
terstellar brightness at 100 /im for each map of our sample. 
The noise level of ea ch map and the CIB fluctuatio n level 
(0.09 MJy sr -1 - see IMiville-Deschenes et al.l (|2002h ) were 
removed quadratically from a. The solid line represents the 
two regimes given by equations [71 and [U 



the power spectrum flattens due to the combination of the 
noise, point sources and the CIB. These three components 
have a power spectrum signature that is flatter than the 
cirrus emission. 

As stated in § [2j the power spectrum of the noise 
is estimated by taking the power spectrum of the dif- 
ference between two HCON images. An example of a 
noise map and its power spe ctrum is shown in Fig. 1 of 
IMiville-Deschenes et al.l ([20021 ) . The power spectrum of the 
noise is well described by a k~ x power spectrum over most 
of the k range. 

To subtract the contribution of bright point sources 
we prefered to removed them directly in the IRIS maps 
prior to compute the power spectrum. To do so we 
used the point source extraction method described by 
IMiville-Deschenes fc Lagachd (|2005| ). For the CIB we as- 
su med a fiat power spect r um a t the level determined 
bv IMiville-Deschenes et all (l2QQ2h (5.8 x 10 3 Jy 2 sr -1 at 
100 /im). Like IMiville-Deschenes et al.l (|2002f ) we made the 
assumption that the effective beam of the IRIS maps is 
Gaussian with a FWHM of 4.3 arcminutes. 



where Po is the power spectrum value at fco = 
0.01 arcmin -1 . The power spectrum of the interstellar emis- 
sion is in general very well fitted by such a power law. 

3.3.1. Spectral index 

The compilation of the power spectrum spectral index mea- 
sured on our sample of 236 maps is shown in Fig. [U The 
most probable spectral index is 7 = —2.9 db 0.2, in accor - 
dance with what was measured bv iGautier et all (|l992f) . 
On the other hand, contrary to IGautier et al.l (|l992f ) who 
did their statistical analysis on only four regions, our larger 
sample allowed us to highlight a significant variation of the 
spectral index with the average brightness of the maps. 
Brighter regions on the sky tend to have a steeper power 
spectrum (see the lower panel of Fig. [3]). The decrease of 
the spectral index with brightness can be approximated by: 

7 = -0.261og 10 (</>)-2.77, (4) 

where < I > is the mean 100 /im brightness of the 12.5° x 
12.5° map. 



3.3. Results 

To obtain the spectral index of the dust emission in each 
IRIS map, we computed the power spectrum of the point 
source subtracted map and the power spectrum of the cor- 
responding noise map. We then subtracted the noise power 
spectrum, divide the result by the Gaussian PSF and re- 
move the CIB contribution. 

Following what was done by IGautier et al.l |l992) we 
have fitted the power spectrum of the Galactic emission in 
the range between k = 0.004 and k = 0.08 arcmin -1 using 
a power law: 

PW = P (£) 7 (3) 



3.3.2. Normalisation 

Regarding the nor malization of the power spectrum 
IGautier et all (|l992f ) estimated that P = 1.4 x 10 6 < I > 3 , 
where < / > is the mean 100 /im brightness at the 12.5° 
scale, given in MJy sr -1 , and Po is in Jy 2 sr -1 . In Fig. [H 
we present the variation of Pp with < / > f or ou r sample, 
together with the relation of IGautier et al.l (|l992l ) (dashed 
li ne). The result of ou r analysis showed that the relation 
of IGautier et al.l (jl992) generally overestimates the fluctu- 
ation level at a given mean brightness. This discrepancy 
could be partly attributed to the fact that we used better 
calibrat ed IRAS data c ompa re to the early IRAS product 
used bv lGautier et al.l (|l992l ). On the other hand it is im- 
portant to point out that the normalization relation given 
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Fig. 6. Left: A typical IRIS map at 100 /im. Center: A classical £Bm map with same power spectrum as the IRIS 
map shown on the left. This fBm has Gaussian brightness fluctuations at all scale and everywhere. Right: A modified 
version of the classical fBm with only positive values and the same average, standard deviation and skewness values as 
the IRIS map. The fBm has been modified in order to reproduce the a (I) oc< I > relation which results in stronger 
brightness fluctuations in brighter regions. 



by iGautier et al.l (|l992l ) is compatible with the fact that 
they used only two faint (< I >~ 4 MJy sr _1 ) and two 

bright regions (< / >~ 30 MJy sr _1 ). 

Using a much larger sample than IGautier et ail (l992) 
we found that the normalization of the power spectrum is 
better described by two regimes (solid line in Fig. on 
- 80% of the sky P(0.01) is proportional to < I > 2 ± 01 
instead of < J > 3 . For < / >< 10 MJy sr" 1 we find 

P(0.01) = 2.7 x 10 6 < I > 2 - 0±0 - 1 . (5) 

and for < / >> 10 MJy sr _1 we find 

P(0.01) = 2.0 x 10 5 < / > 3 - 1±0 - 1 . (6) 

A similar trend was also observed by iJeong et al.l (|2005l ) 
using ISOPHOT observations (see their Fig. 12). 

Another way of looking at the P(0.01)— < / > rela- 
tion is to plot the standard deviation o\ of each map as 
a function of < / > (see Fig. [5]). Here we made sure to 
remove quadratically the contribution of the instrumental 
noise and the CIB to each values of a\. Like for the power 
spectrum normalisation, two regimes are apparent in the o\ 
vs < / > relation with a transition around < / >= 10 MJy 
sr _1 . Separating the data sample in two gives the following 
fits, for < / >< 10 MJy sr" 1 : 

a 2 L = 0.12 < / > 2 - 0±0 - 1 (MJy 2 sr~ 2 ) (7) 
and for < / >> 10 MJy sr _1 

a 2 L = 0.011 < I > 3 - 1±0 - 1 (MJy 2 sr" 2 ) (8) 
3.4. Interpretation 

3.4.1. Spectral index and density structure 

Most spectral indexes measured here fall in the range be- 
tween 7 = — 3.6 to 7 = —2.5 . This is compatible with what 
was found by Wright (1998) (7 = -3) in a power spectrum 
analysis of the DIRBE data at 60, 100, 140 and 240 /im 
on scales greater than 40 arcminutes. Other studies also 



estimated the equivalent of the power spectrum spectral 
index of the 100 urn emission using t he area-perimeter re- 
lation. Following IStutzki et al.l (|l998f ) there is a direct re- 
lation between the fractal dimension D deduced from the 
area-perimeter relation (p ex a D ^ 2 ) and the spectral index 
7 of the power spectrum (P cx k 1 ): 

7 = 2D - 6. (9) 

iBazell fc Desert! (|l988h and iDickman et~aT1 (|l990h mea- 
sured D — 1.25 ± 0.05 (corresponding to 7 = —3.5 ± 0.1) 
for relatively bright regions at 100 /im (~ 10 MJy sr _1 ). 
These values of 7 are rather high (steep power spectrum) 
but not incompatible with our analysis. In a si milar analy- 
sis on fainter regions IVogelaar fc Wakkerl ([l994) measured 
D = 1.45 ± 0.1, corresponding to 7 = —3.1 ± 0.2 which is 
very similar to our results. 

The interpretation of these results in terms of the den- 
sity structure of the ISM and the comparison with what is 
deduced from gas tracers can only be done by considering a 
model for the 100 /im emission. At this wavelength the in- 
terstellar emission is dominated by the grey body emission 
from big d ust grains at therma l equilibrium with the radi- 
ation field ([Desert et all Il990t). For a given model of the 
composition of dust grain ([Desert et al.Ul990l : iLi fc DraineL 
l200lh the conversion of 100 /im brightness to gas column 
density depends on the gas/dust mass ratio (known to be 
rather constant in the ISM) and on the big grain equilib- 
rium temperature. The big grain equilibrium temperature 
is related to the local radiation field strength and spectrum 
which depends on the presence or not of nearby heating 
sources and on the extinction. Variation of the dust equi- 
librium temperature can also occur locally due to variation 
of the grain structure which can affect their emissivit}{3. 

In diffuse regions of the sky at high Galactic lati- 
tudes, far from star- forming regions, clouds are optically 

1 Aggregation of small dust particles on bigger grains can in- 
crease significantly their area / mass ratio which le a ds to a more 
efficient cooling of the grains. See iStepnik et al.l ([2003) for a 
striking example of this effect. 
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thin to stellar radiation, and the radiation field is uniform 
which result in very limited variations of dust equilibrium 
temperature. Localised variations of the dust grain tem- 
perature were observed in cirrus c louds ^ Bernard et all 
1999 ) , but overall several s tudies (|Boulanger fc Peraultl 



19881 : iBoulanger et all Il996l ) showed a strong correlation 



between the 100 jam micron and the hydrogen column den- 
sity which is in favor of a rather uniform gas / dust ratio and 
limited variations of the dust temperature. Based on these 
results we believe that the 100 jam micron can not be used 
as a perfect surrogate for gas column density but overall it 
does not introduce a systematic bias in the determination 
of the column density power spectrum. In this context, and 
considering that the power spectrum of the column density 
gives directly the power s pectru m of the density structure 
([Miville-Deschenes et all l2QQ3bh . we consider that the 7 
values measured here in the diffuse regions (for < / > < 
10 MJy sr _1 ) are typical of the spectral index of the density 
field in three dimensions in the solar neighborhood. 

In bright regions (< I > > 10 MJy sr _1 ) the power 
spectrum is observed to be significantly steeper (Fig [3]). A 
steepeni ng of the power sp ectrum with brightness was also 
found bv lKiss et all (|2003h on 90-200 fim ISOPHOT obser- 
vations. We observed that this steepening coincides with 
a departure from the a oc< I > relation (see Fig. [4j) and 
a systematic increase of the skewness and kurtosis of the 
brightness fluctuations (see Fig. [TT]) . This variation of 7 
with < I > might reflect, at least in part, local variations 
of the density power spectrum. It could also be attributed 
to the effect of gravity or anisotropic radiation fields which 
would both increase the large scale coherence (like in star 
forming regions or at Galactic scale in the plane) and there- 
fore steepen the power spectrum of the 100 /im emission. 
The effect of extinction in these parts of the sky will also 
induce important dust temperature variations which will 
affect the power spectrum. These effects will certainly have 
an impact on the observed brightness fluctuations and make 
the interpretration of the power spectrum difficult. 

To compare the results obtained here with statistical 
analysis of gas emission, it is important to point out that 
the 100 jam emission traces all the interstellar components: 
atomic, molecular and ionized gas. Gas observations (21 
cm and CO) showed that diffuse and dense gas have dif- 
ferent density structure. In general regions of cold and 
dense gas have stronger small scale fluctuations (and there- 
fore a flatter densi ty power spectrum - see for instance 
IStutzki et al.l (|l998l ) who measured 7 = —2.8 on CO obser- 
vations of the Polaris flare, a dense cirrus cloud with a sig- 
nificant fraction of molecular gas ) than regions of diffuse gas 
([Miville-Deschenes et al.l ([2003a] ) measured 7 = —3.6 ± 0.2 
on 21 cm observations of a diffuse cirrus cloud). This be- 
havi or seems also in accordance with numerical simulations 
(see lAudit fc Hennebelld ([2005) for instance). 

This could reflect the fact that molecular tracers like 
CO reveal only the density structure of dense regions and 
miss the more diffuse and large scale structure of molecu- 
lar clouds. In that respect dust emission would be a more 
reliable tracer of the global density structure of the ISM as 
the observed spectral index at 100 jam is a weighted mean 
of the various contributions from dense and diffuse media 
on the line of sight. One striking example of that comes 
from the direct comparison of infrared and 21 cm emission 
in high-latitude clouds where dust emission usually shows 
stronger brightness fluctuations at small scale than Hi, this 



being attributed to th e presence of l ocaliz ed molecular re- 
gions (see for instance I Joncas et al.l (|l992f )). 



3.4.2. The a oc< / > relation 

Here we would like to propose one interpretation for the 
scaling of the power spectrum normalization P(0.01) (or 
equivalently the brightness standard deviation <jl) with av- 
erage brightness < I > (equations [5] and [7]). 

Lets consider a three-dimensional scalar field e, which 
could be a dust emissivity field, of size LxLon the plane 
of the sky and of depth H. In the case of constant dust 
temperature and neglecting opacity effects, the projection 
of this 3D field on 2D would correspond to a dust emission 
map: 



Jo 



I(x,y) = / e(x,y,z)dz. 



(10) 



What would be the relation between the average and stan- 
dard deviation of the projected brightness for such a field ? 
The average brightness of the dust map is simply 



< / >=< e > H 



(11) 



where < e > is the volume averaged of e. 

On the other hand to compute the standard deviation 
of the dust map (07) one should consider that brightness 
fluctuations are added quadratically along the line of sight. 
Each slice dz of the cube contributes o\dz to erf: 



H 



a? dz 



(j 2 c H. 



Combined with Eq. [TT] this leads to 
af oc< I > 



(12) 



(13) 



which is not what is observed. 

Cartesian coordinates used in the previous demonstra- 
tion is in fact not a realistic representation when it comes 
to estimate 07. In fact, with the increase of the depth of 
the line of sight z, a given angular scale 6 on the sky cor- 
responds to increasing physical size I : tan(0) = l/z. 

In the integral of Eq. [12] we made the assumption that 
a e is the standard deviation at scale L. In fact we should 
have taken into account that a e varies with z. For a power 
spectrum following a power law (P(k) ex k 1 ) the vari- 
ance of brightness flu ctuations as a function of scale I is 
([Brunt fc Heverl . 120021 ): 



Gi (X 



1-1/2- 



Therefore we can rewrite Eq. [ 



aj(6)(x / r^- 2 dz= / (tan(6>).?)~ 7 ~ 2 efe 
Jo Jo 

which leads to 

oj(0) oc tan(6»)- 7 - 2 iJ-T- 1 

oc tan(<9)~ 7 " 2 < I >" 7_1 . 



(14) 



(15) 



(16) 



There are two important results to extract from the pre- 
vious equation. First the fact that we observed structure 
of increasing physical size with z at a given angular scale 
does not modify the slope of the power law for < 25°: 
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Fig. 7. Cirrus noise at 100 jam for a lm diffraction lim- 
ited telescope (FHWM=33 arcseconds) as a function of iioo 
brightness. Solid line is our est i mate and dotted line is the 
estimate of iHelou fc Beichmanl (jl99Qh . 



both (Ji(Q) and cr e (Z) have the same spectral index (7 — 2)/2. 
Secondly the dependence of 07 on < I > is now related to 
the spectral index 7. For a typical value of 7 = —3 we ob- 
tain cr| oc< / > 2 which is exactly the observed relation for 
< I >< 10 MJy sr- 1 . 

For larger < I > the departure from a simple aj oc< 
I > 2 can be explained partly because the spectral index 
7 decreases but most of all it is important to notice that 
brightness fluctuations become highly non- Gaussian (high 
skewness and kurtosis - see Fig. [TT]) . In this regime the 
standard deviation becomes significantly affected by non- 
Gaussian fluctuations which could be related to localized 
star forming regions in the image. 



4. Cirrus noise 

In this section we use the properties of the power spectrum 
of 100 jam emission to estimate the level of cirrus confusion 
noise of dust emission. 



4.1. Brightness fluctuations 

To estimate the level of fluctuations as a function of scale 
and brightness it is convenient to use the real-space repre- 
sentation given by Eq. [HJ as opposed to a Fourier one. 
Using the value of 07 measured at I = L = 12.5° (see 
equations [71 and [Hj) we can normalize Eq. [TH and estimate 
the level of brightness fluctuation of interstellar dust (in 
MJy sr -1 ) at scale I (in degree) and wavelength A. For the 
two regimes identified in this study (lower and higher than 
< iioo >=10 MJy sr -1 ) we have the following relations: 



_low 
J l,X 



and 



0.35 < J100 > 



I 



-7/2-1 



aff h = 0.10 < I 100 > 



1.5 J A 
^100 



12.5 C 



-7/2-1 



(17) 



(18) 
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Fig. 8. Contrast of interstellar emission: the standard de- 
viation 0l is the standard deviation of interstellar fluctua- 
tions at the 12° scale. Left: Contrast values as a function 
of average brightness. Right: Histogram of contrast val- 
ues. According to Eq.[7J the contrast is close to constant for 
< I > lower than 10 MJy sr -1 . Outliers at low brightness 
corresponds to regions where there are bright structures on 
a very low level background, like the Magellanic Clouds. 

where < iioo > is in MJy sr -1 , 7 is given by Eq. [H and 
I\/hoo is an estimate of the dust brightness ratio at wave- 
length A and 100 jam. Such rati o can be estimated using a 
grey body model like the one of iDesert et all (jl990h 



A 



h = B u (T d ) 
hoo £100 (T d ) Vl00/xm 



(19) 



where B u (Td) is the Planck function at grain temperature 
Td and (3 is the dust emissivity index. Typical values in 
the diffuse ISM for the dust temperature Td and emissivity 
index /3 are Td = 17.5 K and (5 = 2. In the sub- millimeter 
and millimeter ranges the dust emission departs from this 
simple model and one mig ht want to use a combin ation of 
two dust components (see iFinkbeiner et al.l (|l999f ) for an 
example of such model). 

4.2. Noise per beam 

Equations [T71 and [TBI give the surface brightness fluctuation 
level of dust at any scale, brightness and wavelength. One 
useful specific case to consider is the contribution of cirrus 
noise at the scale of an instrument beam, to estimate the 
effective point sourc e detec tion le vel for instance . Follo wing 
IHelou fc Beichmanl (|l990f ) and iGautier et all ([1992) we 
consider the cirrus noise level at a scale two times the beam 
size b (i.e., I = 26, where b is the beam FWHM). This noise 
level (in mJy/beam) at wavelength A is simply 



lO 9 2M ft 



where Q is the beam in steradian 



V180/ 



(20) 



(21) 



and b is in degrees. The cirrus noise in mJy/beam is then 
given by 

n l r = 3.3 x 10 6 < J100 > 7^ (0.16 6)" 7/2+1 (22) 

^100 
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Fig. 9. Wavelet decomposition of the IRIS map shown in Fig.Oleft. Top: wavelet coefficients map at scales 4, 8, 16 and 
32 pixels (1 pixel=1.5'). Bottom: Histogram of the wavelet coefficients in linear- log. A Gaussian fit to the histogram is 



Gaussian field which is not the case of interstellar emis- 
sion. As it will be described in more details in the next 
section, dust emission has non-Gaussian brightness fluctu- 
ations at all scales. For the sake of data interpretation, con- 
fusion estimate or tests of component separation algorithms 
(including point source extraction) it is useful to produce 
realistic simulations of dust emission maps with proper non- 
Gaussian properties. It is the topic of the next section. 

5. Construction of artificial dust emission maps 

5.1. Fractional Brown ian motion 

Fractional Brownian motions (fBm, also known as Gaussian 
random field s) are often used to simulate interstellar 
emissi on (e.g. JStutzki et alIll998l : lMiville-Deschenes et all 
2003b). By construction such objects can reproduce the 
power spectrum of any image, with the limitation that its 
phase is random. A comparison of a typical IRIS map and 
a fBm with the same power spectrum is given in Fig. [BJ 
The fBm reproduces well the self-similar structure of the 
interstellar emission. On the other hand the fBm seems 
smoother and less contrasted than the observation. 

Because of the Gaussian nature of their fluctuations, 
fBms with positive values only can't have a large contrast. 
If we define the contrast of a map I as Cj = cr(J)/ < I >, 
positive fBms are restricted to C/ < 1/3 as the average 
< I > needs to be greater than 3<r to have only positive 
valuefl This can be compared with the contrast of real 



superposed highlighting the non- Gaussian behavior, 
and 

n h x igh = 1.0 x 10 6 < iioo > 1,5 ^- (0.16 6)- 7/2+1 (23) 

^100 

where 7 is still given by Eq. [H 

These relations can be c ompared with the prescription 
of iHelou fc Beichmanl (|l990f ) often used to estimate cir- 
rus confusion noise. These authors considered the situa- 
tion where the beam of the instrument is given by the 
diffraction limit of a telescope of diameter D. In this case 
b= (1.6A/L>)(180/7r). 

The comparison be tween our estimate of the c irrus noise 
and the one given by IHelou fc Beichmanl (|l990f ) for a lm 
telescope at 100 /xm is shown in Fig [3 The difference is im- 
portant in several aspects. First at low brightness the slopes 
of n(Jioo) are very different, due to the a oc< I > regime 
revealed in our study. The change of slope has the implica- 
tion that our estimate of the cirrus noise in low brightness 
regions is much higher. This is also due to the flattening 
of the power spectrum index 7 at low brightness. On the 
other hand our cirrus noise estimate is much lower than the 
estimate of IHelou fc Beichmanl (|l990f ) at large brightness. 
This effect is caused by the fact that we observe steeper 
power spectra in bright regions. Therefore the brightness 
fluctuations are much smaller when extrapolated to small 
scale. 

4.3. Limitation 



The prescriptions given here to estimate the level of cirrus 
noise is only indicative. It would be strictly correct for a 



2 or equivalently that 99 % of the values are between < I > 
—3a and < I > +3cr. 
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Fig. 10. Wavelet decomposition of the classical fBm map shown in Fig. [6]- center. See Fig. [9] for details. 



IRIS map given in Fig. [HI where we used a(I) = <jl (i-e., 
the standard deviation at a scale of 12.5°). The median 
contrast is 0.3 which indicate the limitation of the use of 
fBms used to simulate realistic infrared dust maps. In ad- 
dition one would notice that the contrast increases signif- 
icantly with brightness and especially for < / > greater 
than 10 MJy sr _1 , in accordance with equations [7] and [8j 
In addition, Cj will increase with scale as a (I) depends on 
scale (according to Eq. [T4]) but not < / > (at high Galactic 
latitude). Based only on the contrast, the use of fBm be- 
comes limited to small (a few degrees) and faint regions of 
the sky. 

5.2. Wavelet decomposition 

Apart from the global contrast limitation, there is a fun- 
damental difference between fBms and observations which 
is related to the non-Gaussian properties of the interstellar 
emission. The brightness fluctuations seen in infrared dust 
maps show contrasted structures, often filaments, that re- 
flect their non- Gaussia n nature . Usin g a wavelet decompo- 
sition of 100 /im maps, iJewelll (|200lf ) showed clearly that, 
contrary to random-phase realisations, the histograms of 
brightness fluctuations at a g iven scale are highly non - 
Gaussian. lAbergel et all ([1996) and lAghanim et all (|2003f ) 
gave also striking examples of that. Wavelet transforms are 
powerful tools to study the statistical moments higher than 
two and estimate the non- Gaussian properties of images 
([Aghanim et all [2003) . They complement the power spec- 
trum analysis which gives only the variation of the second 
moment - the standard deviation - with scale. 

To illustrate that we present in figures [9] and [10] the 
wavelet decomposition obtained using the "a trou" algo- 
rithm (|Starck fc MurtaghL [l998) for the IRIS and fBm 



maps of Fig. [6] As expected the distribution of brightness 
fluctuations at a given scale are Gaussian for fBms (see 
Fig. [TO]) . On the other hand the distribution of wavelet 
coefficients of an IRIS map (see Fig. [9]) follows an asym- 
metrical distribution with exponential wings, which results 
in significant skewness and kurtosis. 

In Fig. [TT] are compiled the skewness and kurtosis of the 
wavelet coefficients found for all 236 maps of our sample, for 
scales from I = 4 pixels (12 arcmin) to 32 pixels (48 arcmin). 
There is a not a strong correlation of the skewness and 
kurtosis with brightness, unlike for the standard deviation 
(Fig. [5]). There is a slight increase of skewness and kurtosis 
at small scale and low brightness that should be attributed 
to noise and CIB 0. On the other hand one would note a 
sharp increase of the skewness and kurtosis for < / >> 
10 MJy sr- 1 . 

Finally, apart from the global difference of the wavelet 
coefficient distribution between fBms and observations, one 
would notice that brightness fluctuations at a given scale in 
the IRIS map are generally greater in bright regions of the 
map, contrary to the fBm where the amplitude of fluctua- 
tions at a given scale is independent of position. This just 
reflects the fact that brightness fluctuations are generally 
stronger in bright region, in accordance with equations [7] 
and [8] By construction fBms do not behave like that; the 
amplitude of fluctuations is uniform and independant of 
local variations of the average brightness. 



3 Both processes being more Gaussian than the interstellar 
emission they have the effect of lowering the skewness and kur- 
tosis mostly. This effect is more important at small scale and 
low brightness where the noise and CIB fluctuations are stronger 
with respect to interstellar fluctuations. 




Fig. 12. Wavelet decomposition of the modified fBm map shown in Fig. fright. See Fig. [9] for details. The histogram 
of the wavelet coefficients of the IRIS map are over plotted (dots) showing how well the modified fBm reproduces the 
observed statistics. 
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Fig. 11. Skewness (top) and Kurtosis (bottom) of the 
wavelet coefficients from 1=4 to 32 pixels for each IRIS map 
of our sample. The solid red line is the median in a bin 20 
values sorted by increasing < I >. 



emission maps that would satisfy a given number of as- 
sumptions: 

1. Positive values only 

2. Power spectrum following a power law : P(k) = Ak 1 

3. Greater brightness fluctuations in bright regions (fol- 
lowing <j oc< / >). 

Several attempts have been made to generate non- 
Gaussian maps, especially to produce non-Gaussian re- 
alizations of the Cosmic Microwave B ackground (CMB) 
(|Vio et all 1200 ll . 120021 : iRocha et all l2005f ) . In general these 
methods involve a modification of a Gaussian realization 
but they don't comply with all our requirements. Especially 
in the context of the CMB there is no need to produce maps 
with only positive values and to have stronger brightness 
fluctuations in bright regions. 

Here is how we proceeded to construct a non- Gaussian 
fBm map with the same statistical properties as an IRIS 
map / and with only positive values. The method pre - 
sented here is similar to the one used by lElmegree n1 (l2002h : 
iBrunt fc Heverl (|2QQ2h . First we generate a classical fBm F 
with same standard deviation as map /. We add an offset 
to force all values to be positives and to be as close as pos- 
sible to the average of /. We then create a modified fBm 
F' such that: 



5.3. Non-Gaussian fBm 

In this section we propose a method to modify fBms such 
that they better reproduce the statistical properties of the 



F' = AF^ +C 



(24) 



and such that F' has the same average, standard deviation 
and skewness than / and only positive values. The parame- 
ters A, ip and C are thus constrained iteratively so to match 
infrared emission. We want to produce simulations of dust the first three statistical moments of the map /. Such a 
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modified fBm is shown in Fig. [6] and its wavelet decompo- 
sition in Fig. [T2J The effect of the exponent ip is to pro- 
duce greater fluctuations in bright regions, in accordance 
with the observations. The wavelet decomposition is quite 
illustrative in that respect. In particular the histograms of 
wavelet coefficients are almost impossible to discriminate 
from the ones of the IRIS map. 

The key parameter in this transformation is ip which 
control the amount of non-Gaussianity introduced in the 
map. We did such a simulation for each of the 236 maps 
of our sample and the results are summarized in Fig. [131 
One would note the scaling of ip with skewness x, which 
clearly indicate the impact of that parameter on the non- 
Gaussianity of the map. The median value of ip is ~ 2 (see 

Fig. us]). 

Even though it reproduces very well several statistical 
properties of the observed emission, it is important to note 
that one important limitation of this method is that it pro- 
duces only isotropic fluctuations and failed to reproduce to 
filamentary structure of the ISM. 

5.4. Extend an image to smaller scales 

To estimate the capabilities and performances of some in- 
struments at observing diffuse interstellar emission it is of- 
ten needed to extrapolate low resolution observations to 
smaller angular scales. One example of that would be the 
estimate of the diffuse emission structure that will be ob- 
served at a scale of 8 arcsec by Herschel-PACS given the 
IRAS data at 5 arcmin resolution. In this context it is use- 
ful to produce constrained realisations of the interstellar 
diffuse emission based on low resolution observations. The 
statistical analysis presented here provides a theoretical ba- 
sis (a ex </>) on which to add small scales using larger 
scale informations. 

Given a low resolution map Io characterized by a beam 
£?o, a higher resolution I\ of beam B\ can be computed 
using the following method: 

h = I + A (B 1 F - B F) it + C (25) 

where F is a classical positive fBm map, and A, ip and C 
are determined as described in the previous section. In the 
previous equation the difference {B\ <g> F — Bo® F) extends 
the structure up to the higher resolution beam B\. The 
modulation of the small scale fluctuations by Jq assures 
that they follow the a oc< / > relation. 

6. Conclusion 

In this paper we presented an analysis of the power spec- 
trum and wavelet decomposition of the IRIS/IRAS 100 /xm 
emission over 55 % of the sky. The main goals of this work 
were 1) to review and extend the study of iGautier et al.l 
(1992) using better calibrated IRAS maps, estimates of the 
noise and CIB contributions and a larger sample, 2) provide 
a more precise prescription for cirrus noise and 3) suggest 
a technique to simulate dust emission map with proper sta- 
tistical properties. 

We found an average spectra l index (< 7 >= — 2.9±0.2) 
compatible with IGautier et al.l (|l992l ) (7 = —3) but with 
a significant variation from 7 = —3.6 to 7 = —2.5. 
Considering that 100 /am emission is a relatively reliable 
tracer of column density in faint regions, these values of 7 
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Fig. 13. Top: Histogram of all the exponent i/j used to 
simulate realistic fBm for each IRIS map of our sample. 
Bottom: Exponent ip as a function of the skewness \ °f 
each IRIS map. The straight line is ip = 2.6 x°" 8 - 



should be representative of the spectral index of density in 
three dimensions in the local interstellar medium. The com- 
parison with other tracers leads to the conclusion that there 
is most probably a significant contribution from cold gas to 
the brightness fluctuations observed. We also found a slight 
variation of 7 with < I > which could be explained by the 
impact of gravity or spatial variations of dust temperature 
in star forming regions. 

We also found that the amplitude of the brightness 
fluctuations were generally overestimated by Gautier et al. 
(1992). In regions with a 100 fim average brightness < I > 
lower than 10 MJy sr _1 the brightness fluctuation level 
is proportional to < / > and not < I > 1,5 as stated by 
IGautier et al.l (p~992). We showed that this behavior can be 
explained by the fact that the brightness fluctuation level 
observed at a given angular size on the sky is the sum of 
fluctuations of increasing amplitude with distance. 

This detailed description of the power spectrum proper- 
ties of the 100 /im emission allowed us to determine a new 
prescription of the cirrus confusion noise in the far-infrared 
and sub-millimeter as a function of column density and 
scale. On the other hand we stressed that this cirrus noise 
estimate relies on the hypothesis of Gaussian fluctuations, 
which is clearly not the case for interstellar emission. In that 
context we proposed a method to modify Gaussian random 
fields such that it reproduces the power spectrum but also 
the level of non-Gaussianity observed which is related to 
the fact that bright regions have stronger brightness fluctu- 
ations than faint ones. Such images could be used to tests 
component separation algorithms (including point source 
extraction) that have to deal with non- Gaussian compo- 
nents. The main limitation of the technique we propose is 
that it does not reproduce the obvious filamentarity seen 
in observations. 



Some of the results in t his paper have been obtained us- 
ing the HEALPix package (jGorski et al.U2005f ). This work 
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